## This program generates Figure 2 from "Financial Constraints, Sectoral Heterogeneity,
##  and the Cyclicality of Investment" by Cooper Howes

##### Step 0: Preliminaries #####

## Clear variables and data
rm(list=ls())

## Load packages
library(sandwich)
library(tidyverse)
library(glue)
library(dynlm)
library(seasonal)

## Set maximum length of response horizon in quarters (default:16)
numres <- 16 

## Set number of standard deviations to show in errors (default: 1.65, i.e. 90% confidence intervals)
se_scale <- 1.65

## Set scale of shock in percentage points (default: 1)
scale <- 1 

## Set number of lags to use in regressions
lags <- 8

## Specify shocks to use in regression (default: "rr", or Romer & Romer (2004)-style narrative shocks)
shocktype <- "rr"

## Regression includes linear time trend, seasonal term, lags of DV, and lags of shock
f <- formula( glue(paste("L(y,-i) ~ trend(y) + season(y) + {shocktype}_shock", 
                        "+ L({shocktype}_shock,1:lags)  ", 
                        " +L(y,1:lags) " )))


## Set beginning and ending dates for regressions in (yyyy,q) format
begin <- c(1970,1) ## Default: c(1970,1)
fin   <- c(2008,4) ## Default: c(2008,4)
                

## Load data from 1966-2021
data <- ts(read.csv("Aggregate_data.csv"),start=c(1966,1),end=c(2021,4),frequency=4)

## Specify which variables to plot
varnames <- c("rsales_all","rsales_dur","rsales_nondur",
              "rnppe_all","rnppe_dur","rnppe_nondur")

outcomes <- data[,varnames]

plotvars <- 1:length(varnames)

## Create empty arrays to hold coefficients and standard errors
gamma    <- matrix(0,nrow=(numres+1),ncol=length(varnames))
gamma_se <- matrix(0,nrow=(numres+1),ncol=length(varnames))


##### Step 1: Create top panels (IRFs for sales and capital stocks) #####


## Outer loop over IRFs to be estimated
for(j in plotvars){
  
  n <- varnames[j]
  
  y <- log(data[,n])
  

  ## Inner loop over each quarter in the response horizon
  for(i in 0:numres){
    
    ## Need to calculate date of last shock used so all data end in 2008. For
    ##  example, if the data end in 2008, then the latest shock used when estimating
    ##  the response 13 quarters after the shock will need to be in the third
    ##  quarter of 2005, since (2005Q3 + 13 quarters = 2008Q4).
    lastshock_y <- fin[1]-floor(i/4)
    lastshock_q <- 4-i%%4
    
    ## Estimate regression and extract Newey-West standard errors
    eq <- dynlm(formula=f,start=begin,end=c(lastshock_y,lastshock_q),data=data)
    newey <- NeweyWest(eq,prewhite=FALSE)
    
    ## Save coefficients and standard errors
    gamma[(1+i),j] <- coef(summary(eq))[paste0(shocktype,"_shock"),"Estimate"]
    gamma_se[(1+i),j] <- sqrt(newey[paste0(shocktype,"_shock"),paste0(shocktype,"_shock")])
    
  }
  
  
}

##### Step 2: Create bottom panels (IRFs for sales/capital stock "gaps") #####


## Create empty arrays to hold coefficients and standard errors
gapmat <- matrix(0,nrow=(numres+1),ncol=2)
gapmat_se <- matrix(0,nrow=(numres+1),ncol=2)

## Calculate sales and NPPE "gaps" 
salesgap <- log(data[,"rsales_dur"]) - log(data[,"rsales_nondur"])
nppegap <- log(data[,"rnppe_dur"]) - log(data[,"rnppe_nondur"])

## Loop over each quarter of the response horizon
for(i in 0:numres){
  
  ## Need to calculate date of last shock used so all data end in 2008. For
  ##  example, if the data end in 2008, then the latest shock used when estimating
  ##  the response 13 quarters after the shock will need to be in the third
  ##  quarter of 2005, since (2005Q3 + 13 quarters = 2008Q4).
  lastshock_y <- fin[1]-floor(i/4)
  lastshock_q <- 4-i%%4
  
  ## Estimate response of sales gap to MP shock and save coefficients and SEs
  y <- salesgap
  eq <- dynlm(formula=f,start=begin,end=c(lastshock_y,lastshock_q),data=data)
  newey <- NeweyWest(eq,prewhite=FALSE)
  
  gapmat[(1+i),1] <- coef(summary(eq))[paste0(shocktype,"_shock"),"Estimate"]
  gapmat_se[(1+i),1] <- sqrt(newey[paste0(shocktype,"_shock"),paste0(shocktype,"_shock")])
  
  
  ## Estimate response of NPPE gap to MP shock and save coefficients and SEs
  y <- nppegap
  eq <- dynlm(formula=f,start=begin,end=c(lastshock_y,lastshock_q),data=data)
  newey <- NeweyWest(eq,prewhite=FALSE)
  
  gapmat[(1+i),2] <- coef(summary(eq))[paste0(shocktype,"_shock"),"Estimate"]
  gapmat_se[(1+i),2] <- sqrt(newey[paste0(shocktype,"_shock"),paste0(shocktype,"_shock")])
  
  
}


#####  Step 3: Plot top panels (sales and capital stock IRFs) #####


qtrs <- c(0:numres) ## Generates x-axis sequence for plots


## Set parameters for plot windows
w       <- 2*length(varnames)         ## Window width
h       <- length(varnames)         ## Window height

## Specify which IRFs will be plotted
varnames2 <- c("rsales_dur","rsales_nondur","rsales_all","rnppe_dur","rnppe_nondur","rnppe_all")


## Create new window
windows(width=12,height=5)
par(oma=c(4,2,0,0))
m <- matrix(1:2,nrow=1,byrow = TRUE)
layout(mat = m)

## Plot sales IRFs

# Durable sales + SEs
g1 <- plotvars[2]
irf1<-scale*gamma[,g1]
irf2<-scale*gamma[,g1]+se_scale*gamma_se[,g1]
irf3<-scale*gamma[,g1]-se_scale*gamma_se[,g1]

# Nondurable sales + SEs
g2 <- plotvars[3]
irf4<-scale*gamma[,g2]
irf5<-scale*gamma[,g2]+se_scale*gamma_se[,g2]
irf6<-scale*gamma[,g2]-se_scale*gamma_se[,g2]

# Total sales 
g3 <- plotvars[1]
irf7 <- scale*gamma[,g3]

## Automatically generate axis ranges
y_range1 <- range(irf1,irf2,irf3,irf4,irf5,irf6,irf7)

## Plot durable sales coefficient estimates with filled SE bands
par(mar = c(2,2,2,2))
plot(qtrs,irf2,type="l",lty=2,col="blue",lwd=2,ylim=y_range1,ann=FALSE,cex.main=1.25,axes=F,frame=T)
lines(qtrs,irf3,type="l",lty=2,col="blue",lwd=2)
polygon(c(qtrs,rev(qtrs)), c(irf3,rev(irf2)), col=rgb(red=0,blue=1,green=0,alpha=0.5),border=NA)
lines(qtrs,irf1,type="b",pch=1,lwd=3,col="blue")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5,at=pretty(y_range1),lab=paste0(pretty(y_range1)*100, "%"),las=TRUE)

## Add nondurable coefficients and SEs
lines(qtrs,irf5,type="l",lty=2,col="red",lwd=2)
lines(qtrs,irf6,type="l",lty=2,col="red",lwd=2)
polygon(c(qtrs,rev(qtrs)), c(irf6,rev(irf5)), col=rgb(red=1,blue=0,green=0,alpha=0.5),border=NA)
lines(x=qtrs,y=irf4,type="b",col="red",lwd=3,pch=2)

## Plot aggregate sales and add a line at y=0 and a title
lines(x=qtrs,y=irf7,type="b",col="orange",lwd=3,pch=5)
abline(a=0,b=0,col="black",lwd=1)
title(main="Real sales", col.main="blue", font.main=2,cex.main=1.25)


## Plot NPPE IRFs

## Durable NPPE + SEs
g4 <- plotvars[5]
irf8<-scale*gamma[,g4]
irf9<-scale*gamma[,g4]+se_scale*gamma_se[,g4]
irf10<-scale*gamma[,g4]-se_scale*gamma_se[,g4]

## Nondurable NPPE + SEs
g5 <- plotvars[6]
irf11<-scale*gamma[,g5]
irf12<-scale*gamma[,g5]+se_scale*gamma_se[,g5]
irf13<-scale*gamma[,g5]-se_scale*gamma_se[,g5]

## Total NPPE
g6 <- plotvars[4]
irf14 <- scale*gamma[,g6]

## Automatically generate axis ranges
y_range2 <- range(irf8,irf9,irf10,irf11,irf12,irf13,irf14)


## Plot durable NPPE coefficient estimates with filled SE bands
par(mar = c(2,2,2,2))
plot(qtrs,irf9,type="l",lty=2,col="blue",lwd=2,ylim=y_range2,ann=FALSE,cex.main=1.25,axes=F,frame=T)
lines(qtrs,irf10,type="l",lty=2,col="blue",lwd=2)
polygon(c(qtrs,rev(qtrs)), c(irf10,rev(irf9)), col=rgb(red=0,blue=1,green=0,alpha=0.5),border=NA)
lines(qtrs,irf8,type="b",pch=1,lwd=3,col="blue")
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5,at=pretty(y_range2),lab=paste0(pretty(y_range2)*100, "%"),las=TRUE)

## Add nondurable coefficients + SEs
lines(qtrs,irf12,type="l",lty=2,col="red",lwd=2)
lines(qtrs,irf13,type="l",lty=2,col="red",lwd=2)
polygon(c(qtrs,rev(qtrs)), c(irf13,rev(irf12)), col=rgb(red=1,blue=0,green=0,alpha=0.5),border=NA)
lines(x=qtrs,y=irf11,type="b",col="red",lwd=3,pch=2)

## Plot aggregate NPPE and add a line at y=0 and a title
lines(x=qtrs,y=irf14,type="b",col="orange",lwd=3,pch=5)
abline(a=0,b=0,col="black",lwd=1)
title(main="Real capital stock", col.main="blue", font.main=2,cex.main=1.25)

## Add legend at the bottom
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend(x="bottom",inset=0,legend=c("Durable","","Nondurable","","All"),
       lwd=2,pch=c(1,NA,2,NA,5),col=c("blue",NA,"red",NA,"orange"),x.intersp=1,cex=1.25,horiz=TRUE)


## Save empirical IRFs for comparison with the model in Figure 6
qfr_irfs <- cbind(irf8,irf11,irf14)
colnames(qfr_irfs) <- c("dur","nondur","total")
write.csv(qfr_irfs,"irfs/qfr_irfs.csv")


#####  Step 4: Plot bottom panels (sales and capital stock gap IRFs) #####

## Create new window
windows(width=12,height=5)
par(oma=c(4,2,0,0))
m <- matrix(1:2,nrow=1,byrow = TRUE)
layout(mat = m)

## Sales gap IRFs
irf1<-scale*gapmat[,1]
irf2<-scale*gapmat[,1]+se_scale*gapmat_se[,1]
irf3<-scale*gapmat[,1]-se_scale*gapmat_se[,1]

## NPPE gap IRFs
irf4<-scale*gapmat[,2]
irf5<-scale*gapmat[,2]+se_scale*gapmat_se[,2]
irf6<-scale*gapmat[,2]-se_scale*gapmat_se[,2]

## Automatically generate axis ranges
y_range1 <- range(irf1,irf2,irf3)
y_range2 <- range(irf4,irf5,irf6)

## Create plot for sales gap IRFs
par(mar = c(2,2,2,2))
plot(qtrs,irf1,type="l",lty=1,col="blue",lwd=3,ylim=y_range1,ann=FALSE,cex.main=1.25,axes=F,frame=T)
lines(qtrs,irf2,type="l",lty=2,col="red",lwd=2)
lines(qtrs,irf3,type="l",lty=2,col="red",lwd=2)
abline(a=0,b=0,col="black",lwd=1)
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5,at=pretty(y_range1),lab=paste0(pretty(y_range1)*100, "%"),las=TRUE)
title(main="Real sales gap", col.main="blue", font.main=2,cex.main=1.25)

# Create plot for NPPE IRFs
plot(qtrs,irf4,type="l",lty=1,col="blue",lwd=3,ylim=y_range2,ann=FALSE,cex.main=1.25,axes=F,frame=T)
lines(qtrs,irf5,type="l",lty=2,col="red",lwd=2)
lines(qtrs,irf6,type="l",lty=2,col="red",lwd=2)
abline(a=0,b=0,col="black",lwd=1)
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5,at=pretty(y_range2),lab=paste0(pretty(y_range2)*100, "%"),las=TRUE)
title(main="Real capital stock gap", col.main="blue", font.main=2,cex.main=1.25)

## Add legend for bottom panels
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend(x="bottom",inset=0,lty=c(1,2),legend=c("Coefficient Estimate","90% Confidence Interval"),
       lwd=2,col=c("blue","red"),cex=1.25,horiz=TRUE)

